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Levy Processes For Finance: 
An Introduction in R 


D.J. Manuge 


Abstract 

This brief manuscript provides an introduction to Levy processes and their applications in finance as the 
random process that drives asset models. Characteristic functions and random variable generators of popular 
Levy processes are presented in R. 


1. Introduction 

Over the past few decades analysis has shown that market data is inconsistent with some of the underlying 
assumptions of the Black-Scholes model. For example, if asset price changes remain large while time periods 
shrink, one cannot assume that prices are continuous. Under this observation, Cox and Ross assume prices 
follow a pure jump process so that at any timestep the asset price results in a positive jump or negative 
drift [7]. Expanding on this theory, Merton introduced jumps superimposed on a continuous price process, 
whose parameters can be chosen to account for fat-tail returns observed in real markets [22]. Additionally, 
practitioners require models with accurate implied volatility surfaces for risk management. While the 
surface can be explained without a jump model, the smile becomes augmented for short maturity options, 
which is better described by the presence of jumps ED.Q Moreover, distributions of returns exhibit skewness 
and leptokurtosis. Clearly, as more information about market structure is understood, improvements to 
parametric forms of asset prices are neccessary. For the class of finite jump models, Kou extends Merton’s 
model by allowing jump sizes to follow an asymmetric double exponential distribution |16j . The model 
accounts for both the volatility smile and asymmetric leptokurtic returns. Most financial models fall into 
the class of infinite jump models, where many developments have been made. In 1987, Madan and Seneta 
suggest that increments of log-prices follow a (symmetric) variance gamma (VG) distribution, while providing 
statistical evidence using Australian stock market data [niEi]. a The VG distribution is a special case of 
the generalized hyperbolic (GH) distribution, which was initially proposed by Barndorff-Nielson to model 
the grain-size distribution of sand as it travels from a source to a deposit [5|. Other cases of the GH 
distribution have been considered for asset price modelling; Eberlein and Keller suggest that increments of 
log-prices follow a hyperbolic distribution [5|. Alternatively, Barndorff-Nielson propose the normal inverse 
Gaussian (NIG) distribution for log-prices in |3|. Finally, Eberlein and Prause extend their focus to the entire 
class of GH distributions, providing a thorough statistical treatment as well as applications to stochastic 
volatility |SJ[5^. Carr, Geman, Madan, and Yor proposed the (generalized) tempered stable process. |5]. 
Coined the CGMY process, it coincides with the Levy measure of a non-normal a-stable Levy process 
{a € (0,2)) multiplied by an exponential factor. Another infinite jump model, which arose from the theory 
of orthogonal polynomials, is given by the Meixner process [29l[30|. All of the processes mentioned thus far 
share a common trait; they are Levy processes. It is obvious that by modelling the asset dynamics with a 
general Levy process, a large range of special cases are accessible for consideration. Not only does this 
general class of process sufficiently mimic the implied volatility surface of real data, but as mentioned they 


^Known as the smile effect of volatility, due to its graphical portrayal of happiness. 

^In non-jump models, as T —> 0, \n{St/St—i) ~ Normal(/i, tr), when in reality this is not the case. 

® Later extended to the asymmetric variance gamma model in m 

^Which include (but are not limited to) the financial cases mentioned above. 
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can be parametrized to exhibit skewness, kurtosis, an absence of autocorrelation in price increments, finite 
variance, aggregational normality, and have an ability to change discontinuously [25l[24l[2^ . 

2. Levy Processes 

Definition 1 (Stochastic Process) A stochastic process X on a probability space (fi, P) is a colleetion 
of random variables (-^t)o<i<oo- 

If Xt G iFt, the process X is said to be adapted to the filtration X, or equivalently. Ft — measurable (1271). 
Stochastic processes have been studied in economic, financial, actuarial, physical, biological, and chemical 
applications m. m. E], m, iia). In finance, their primary use is to represent the price component 
of some asset. Arguably, the most well-known stochastic process is Brownian motion, which describes the 
movement of particles suspended in a fluid ([in],[n]). 

Definition 2 (Brownian Motion) Standard Brownian motion W = (Wt)o<t<oo satisfies the following 
three properties: 

(i) Wo = 0 

(a) W has independent increments: Wt — Wg is independent of Fg, 0 < s < t < oo 
(Hi) Wt — Wg is a Gaussian random variable: Wt — Wg ^ N{0, t — s)V0<s<t<oo 

Property (ii) implies the Markov property (i.e. conditional probability distribution of future states depend 
only on the present state). Property (Hi) indicates that knowing the distribution of Wt for t < r provides 
no predictive information about the process when t > t. Another well-known stochastic process is l|13j'l. 

Definition 3 (Poisson Process) A Poisson process N = (A^t)o<t<oo satisfies the following three proper¬ 
ties: 

(i) No = 0 

(ii) N has independent increments: Nt — Ng is independent of Fg, 0 < s < t < oo 
(Hi) N has stationary increments: P{Nt — Ng < x) = P{Nt-g < x), 0 < s < t < oo 

In reality, stochastic differential equations formulated with only Brownian motion or the Poisson process 
may be ineffective in describing the complex dynamics of an evolving system. Although their path structures 
appear different in simulation, note the similarity in their definitions. By combining their common properties, 
a general process can be established. 

Definition 4 (Levy Process) Let L be a stochastic process. Then Lt is a Levy process if the following 
conditions are satisfied: 

(^) Lo = 0 

(ii) L has independent increments: Lt — Lg is independent of Fg, 0 < s < t < oo 
(Hi) L has stationary increments: P(Lt — Lg < x) = ¥{Lt-g < x), 0 < s < t < oo 
(Hi) Lt is continuous in probability: limt_>s Lj = Lg 

Remark 1 Keller shows that condition (Hi) follows from (i) and (ii), and thus can be omitted (JJ^)- 

A random variable L is said to have an infinitely divisible distribution if there exists a sequence of i.i.d 
random variables such that L = -|- ... -|- Ln^n for each n G It follows that if L is a Levy process, it 

has an infinitely divisible distribution for each t > 0 (irzi)^0 An expression exists which fully characterizes 
the properties stated, giving a parametric structure to this process’ seemingly qualitative definition. 


^This is merely one of many results attributed to Paul Levy (1886-1971) who produced foundational results in the calculus 
of probabilities; at a time when no mathematical theory of probability was available. His doctoral advisors were Jacques 
Hadamard and Vito Volterra. 

®The converse holds as well; all infinitely divisible distributions can be represented as a Levy process. 


2 



Theorem 1 (Levy-Khintchine formula) Suppose that /x S M, cr > 0, and v is a measure concentrated 
on M/{0} such that J^imn{l,x'^)v(dx) < oo. A probability law r] of a real-valued random variable L has 
characteristic exponent ^'(m) := —^ logE[e“'^*] given by, 


= / e“^r?(dx) = for uGR, 

Jr 

iff there exists a triple {p, a, v) such that, 


'P('u) = ipLu + + / (1 — + zua:lL(|3,|<i))jx(cix) 

2 Jr 


( 2 . 1 ) 


( 2 . 2 ) 


for every m € K. 

Furthermore, it follows from the above theorem that there exists a probability space where L = + 

is standard Brownian motion with drift, is a compound Poisson process, and is 

a square integrable martingale with countable number of jumps of magnitude less than 1 (almost surely). 
The characteristic exponent of the process is the same as in the Levy-Khintchine formula. This result is 
more formally known as the Levy-Ito decomposition (EZl)- In canonical form. 


Lt = pitaWt + 


xr]^{ds, dx) 


x{rj^ — i’^){ds, dx). 


0 |x|>l 


0 |a:|<l 


(2.3) 


3. Simulating Assets with Levy Processes 

The path structure of any Levy process is uniquely defined by its triplet {pi,a,v). Many popular asset 
pricing models and their associated driving processes have known triplets m, ini)- This section discusses 
examples of Levy processes that are used in finance to model assets. In the cases which follow we assume 
the price process takes the form. 


S = exp(Lt) 

where Lt is a Levy process. In this paper, we do not discuss details regarding the risk-neutrality of the 
asset. However, we note that the asset must be mean-corrected for it to be risk-neutral. One way this can 
be achieved is by multiplying the above expression by an exponential damping factor which can be solved 
by knowing the distribution of the asset (i.e. the choice Levy process) ([H]). In the code that presents 
the characteristic functions of Levy processes below, risk-neutrality is applied. The ” fAsianOptions” and 
’’Bessel” R packages should be installed and loaded prior to running the scripts in this section. 

Brownian Motion: Suppose we have mean pi and variance with Li ~ Normal(^, cr^), then the triplet 
is given by (^, cr^,0). The density t^Li is given by, 

= (3.1) 

V27rcr 

The canonical decomposition is given by, 

Lt = pit -\- crWt (3-2) 

which has a characteristic exponent of the form, 

'I'(u) = ~ iup. (3.3) 

which can be obtained via: 




martingale is a process L whose conditional expected value at is equal to Ln. 
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BS_CF=function(u,sigma.r,time){ 
drift=r-0.5*sigma~2 

phi=exp(li*drift*time*u-0.5*sigma*sigma*u*u*time) 
return(phi) 

} 

Random variables can be obtained to simulate the path of assets by calling: 

BM=function(mu,sigma, T, N) { 

h=T/N 

t=(0:T)/N 

X=rep(0, N+1) 

X[1]=0 

for(i in 1:N) { X[i+l]=X[i] +mu*h+sigma*sqrt(h)*rnorm(l)} 
return(X) 

} 


Asset price process driven by Brownian motion 



Time 


Figure 1: Sample paths of an asset driven by Brownian motion with = 0.5 and a = 0.5. 

Log-returns are normally distributed in the BSM model; the log-price process follows a standard Brownian 
motion with drift m)- Setting ^ — 0 eliminates the drift component. 

Poisson Process: Let A > 0 and 5i be the Dirac measure for x = 1. Suppose Li ^ Poisson(A), then the 
triplet is given by (0,0, A^i). The density rjLi is given by, 

= Ae“^“. (3.4) 

The decomposition is the cumulative sum of jumps up to time t, which has a characteristic exponent of the 
form, 

4-(u) = A(l-e“). (3.5) 

Random variables can be obtained to simulate the path of assets by calling: 


®5i(T) = 1 when 1 £ A, otherwise (5i(T) = 0. 
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#Poisson generator 
PPgen=function(lambda){ 

X=0 
Sum=0 
flag=0 

while (flag==0){ 

E=-log(runif(1)) 

Sum=Sum+E 

if (Sum < lambda) ■[ X=X+1} else { flag=l} 

} 

return(X) 

} 

#Poisson process 

PP=function(lcmibda, N, T){ 

h=T/N 

t=(0:T)/N 

X=rep(0, N+1) 

l=rep(0,N) 

X[1]=0 

ford in 1:N) { 

1 [i]=PPgen(h*lambda) 

X[i+l]=X[i] + l[i]} 
return(X) 

} 


Poisson process 



Figure 2: Sample paths of a Poisson process with A = 100. 


Compound Poisson Process: Let A > 0 be the intensity of a Poisson random variable N and let 
be random variables with law F for * > 1 (independent of N). Suppose Li follows a compound Poisson 
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distribution. Then the triplet is given by (—A /g<| 2 ,|<i xF{dx), 0, \F{dx)). The density 77 ^^ is not analytically 
tractable. The canonical decomposition is given by, 

Nt 

Lt = Y.^- (3-6) 

i=l 


which has a characteristic exponent of the form, 

■qj{u) = A /(I - F^^)F{dx). (3.7) 


Random variables can be obtained to simulate the path of assets by calling: 

CPP=function(lambda, T, N) { 

h=T/N 

t=(0:T)/N 

X=rep(0, N+1) 

F=rep(0, N+1) 

I=rep(0,N) 

X[1]=0 

ford in 1:N) { 

I [i]=PPgen(h*lambda) 

if (I[i]==0){F [i]=0} else {F[i]=rnorm(l)} 

X[i+l]=X[i] + F[i]} 
return(X) 

} 


Compound Poisson process 



Figure 3: Sample paths of a compound Poisson process with A = 100 and a standard normal jump distribution. 


^Note that the function PPgen is required from before. 
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Merton model: Suppose that log-returns Li are determined by the Merton model; the sum of standard 
Brownian motion (with drift) and a compound Poisson process whose jumps are normally distributed ([22]). 
Then ~ Normal(cr^) for i > 1 and the triplet is (/i,cr^,A x ryj). As expected, the density rjLi is not 
analytically tractable. The density is given by, 



(3.8) 


The canonical decomposition is given by, 


Lt — ^it + uWt + 


(3.9) 




which has a characteristic exponent of the form, 


^(U) = ^ 


(3.10) 


which can be obtained via: 

Merton_CF=function(u,sigma,a,b,lambda, time){ 

jump=lambda*time*(-a*u*li+(exp(u*li*log(l+a)+0.5*b*b*u*li*(u*li-l))-l)) 
phi=BS_CF(u,sigma,r,time)+exp(jump) 
return(phi) 

} 

Random variables can be obtained to simulate the path of assets by calling: C3 

Merton=function(mu,sigma,lambda,mu_xi,sigma_xi, N, T) { 

h=T/N 

t=(0:T)/N 

X=rep(0, N+1) 

F=rep(0, N+1) 
l=rep(0,N) 

X[1]=0 

ford in 1:N) { 

1 [i]=PPgen(h+lambda) 

if (1 [i]==0){F [i]=0} else {F [i]=mu_xi*l [i]+ sqrt(sigma_xi)*sqrt(1[i])*rnorm(l)} 
X[i+l]=X[i] + mu+h+sigma*sqrt(h)*rnorm(l)+F [i]} 
return(X) 

} 


^'^Note that the function PPgen is required from before. 
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Asset price driven by Merton process 



Time 


Figure 4: Sample paths of an asset driven by a Merton process with fi = 0.5, a = 0.75, and A = 1.5. 


Kou model: Suppose that log-returns Li are determined by the Kou model; the sum of standard Brownian 
motion (with drift) and a compound Poisson process whose jumps are double-exponentially distributed ('|16j'). 
Then ~ DExp(p, 0i, @ 2 ) for t > 1 and the triplet is x rj^). Once again, the density rjL^ is not 

analytically tractable. The density is given by, 


=pOie + (1 -p)6'2e®""'l{x>o}■ 

The canonical decomposition is given by, 

Nt 

Lt = pt + aWt + 

i=l 


which has a characteristic exponent of the form. 


'I'(u) 


^ 2 2 
-a u 


iufjL + A 



{l-p)d2 
O 2 + iu 


pOi \ 
9i — iu J 


(3.11) 


(3.12) 


(3.13) 


which can be obtained via: 

Kou_CF=function(u,sigma, lambda, p, theta_l, theta_2, r, time){ 

drift=r-0.5*sigma~2 -lambda*(((p*theta_l)/(theta_2+l))+((1-p)*theta_2/(theta_2+l))-1) 
phi=exp(li*u*drift*time-0.5*sigma"2*u*u*time+time*lambda*((p*theta_l)/(theta_l+u*li)+((1-p)*theta_2)/ 
return(phi) 

} 

Random variables can be obtained to simulate the path of assets by calling: 
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Kou=functioii(mu,sigma,lcmibda,p,theta_l,theta_2, T, N) { 

h=T/N 

t=(0:T)/N 

X=rep(0, N+1) 

F=rep(0, N+1) 

I=rep(0,N) 

X[1]=0 

ford in 1:N) { 

I [i] =PPgen(li*lambda) 
if (I[i]==0){F[i]=0} 
else { 

K=rbinom(l,I[i],p) 

Rl=rgaiimia(l ,K*theta_l ,theta_2) 

R2=rgaiiima(l, (I [i] -K) *theta_l, theta_2) 

F[i]=Rl-R2} 

X[i+l]=X[i] + mu+h+sigma*sqrt(h)*rnorm(l)+F [i]} 
return(X) 

} 


Asset price driven by Kou process 



Figure 5: Sample paths of an asset driven by a Kou process with fi = 0.5, cr = 0.75, A = 1.5, p = 0.5, 9i = 0.25, and 02 = 0.25. 

Jump-diffusion models are generated by processes with a non-zero diffusion component and finite jumps. 
The above examples are all of this type since cr > 0 and v{dx) < oo. If instead cr = 0 and v{dy) = oo, we 
obtain the class of infinite activity pure jump models. The majority of asset pricing models proposed in 
recent years are of this type. 

Variance Gamma Process: Suppose that log-returns Li are determined by a variance gamma price 
process m, m)- The variance gamma process can be obtained by evaluating standard Brownian motion 
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with drift, subordinated by a gamma process in the time variable. Let fj,,a,v > 0 

1 


C= - > 0 

V 


- i/ 
1 

2 ^ 


-1 


— ii2q,2 _i_ —/x2^ 


V 2 

By defining 7 = ~ 1) ~ G{e~^ — 1)] and, 


the triplet becomes ( 7 , 0 , 2 ^'^'^). The density rjL^ is not analytically tractable (page 4, rathgeber, stabler, 
stockl). The canonical decomposition is given by, 



(3.14) 

> 0 

(3.15) 

> 0 . 

(3.16) 

G^t^yodx 

(3.17) 


Lt = tE[Li] + / / x{r]^ — i'^^){ds,dx) 


(3.18) 


which has a characteristic exponent of the form. 


4>(m) = i In [1 — iu^v + \o'‘^vu^\. 

V 2 


(3.19) 


which can be obtained via: 


VG_CF=function(u,sigma, theta, nu, r, time){ 
drift=r+log(l-theta*nu -0.5*sigma*sigma*nu)/nu 

phi=exp(li*drift*time*u) *((l-li*nu*theta*u+0.5*nu*sigma*sigma*u*u)~(-time/nu)) 
return(phi) 

} 


Alternatively, (3.18) can be expressed as the difference between two gamma processes i.e. Lt = G\ — G^. 
The density of a gamma process is given by, 


= 


r(a) 




(3.20) 


where a, 6 , x > 0. By setting a = G and 6 = M for G\ and a = G and & = G for , one can simulate the 
VG process in this way. 

VG=function(sigma, nu, mu, T, N) { 

a=l/nu 

b=l/nu 

h=T/N 

t=(0:T)/N 

X=rep(0, N+1) 

I=rep(0,N) 

X[1]=0 

ford in 1:N) { 

I[i]=rgamma(1,a*h,b) 

X[i+l]=X[i] + mu*I[i]+sigma*sqrt(I[i])*rnorm(l) 

} 

return(X) 

} 


^^Note that the function PPgen is required from before. 
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Asset price driven by VG process 



Time 


Figure 6: Sample paths of an asset driven by a VG process with cr = 0.75, v = 0.5, and fi = 0.1. 


CGMY Process: Suppose that log-returns Li ~ CGMY(C, G, M, Y) where G, G, M are adopted from 
(3.14)-(3.16) and Y <2 ([5]). This process has infinite activity iff Y € [0, 2). The first moment does not 
have an explicit form. Let, 


(j^-Mx (j^-G\x\ 

V (dx) — -f (3.21) 

then the triplet is given by (E[Li], 0, v"'"^{dx)). The density is not analytically tractable. The canonical 
decomposition is given by, 


(3.22) 


(3.23) 

which can be obtained via: 

CGMY_CF=function(u,T,r,C,G,M,Y){ 
omega=-C*cgainma(-Y)*((M-1)“Y-M~Y+(G+li*u)~Y-G~Y) 
tmp=C*T*cgaiiima(-Y) * ( (M-li*u) "Y-M"Y+ (G+li*u) "Y-G"Y) 
phi=exp(1i*u*((r+omega)*T)+tmp) 
return(phi) 

} 


Lt=t'E\Li]+ I I x{ri^ — v'^^){ds,dx) 


which has a characteristic exponent of the form, 

4-(m) = Gr(-Y)[G^ + M^ -{M- iuf -{G + iu)^]. 


^^With a change of variables, this is the tempered stable process. 
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Similar to the VG process, the CGMY process can be simulated using Brownian subordination m- Let, 


G-M 



(3.24) 


Then the process defined by Lt = AYt + WiYt) is equivalent to (3.22[) where Y* is the process with density. 



(3.25) 


and for F > 0, 



(3.26) 


Random variables can be obtained to simulate the path of assets by calling: 

CGMY=function(C,G,M,Y, T,N) { 

h=T/N 

t=(0:T)/N 

A=(G-M)/2 

B=(G+M)/2 

f=function(y){exp(-(B~2-A~2)*y/2)*cgainma(Y)*CGMY_f(y,Y,B)/(cgcmima(Y/2)*2"(Y/2-1))} 
X=rep(0, N+1) 

I=rep(0,N) 

probJ=rep(0,N) 

J=(0:(N-l))*h + h 

for (j in 1:N) {probJ [ [j] ] <-f (J [ [j]] )} 
randf=sample(J,N,probJ,replace=TRUE) 

X[1]=0 

ford in 1:N) { 

I[i]=randf [[i]] 

X[i+l]=X[i] + A*I [i]+sqrt (I [i] ) *rnorni(l) 

} 

return(X) 

} 
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Asset price driven by CGMY process 



Figure 7: Sample paths of an asset driven by a CGMY process with C = 5, G = 25, M = 25, and Y = 1. 


Generalized Hyperbolic Process: Suppose that log-returns Li ~ GH(a,/?, <5, A) ([9], [26], |2]). The 
parameter a > 0 determines the shape, 0 < |/3| < a determines the skewness, /i € K gives the location, 
(5 > 0 is the scaling parameter, and A S K determines the heaviness of the tails. The first moment is given 

by, 

PS^K x+iiO 


E[Li] = n 


where C = S\/a'^ — By defining 

-"■"■Si/ 

\ n 


CKxiO 


-dy + Xe 


('^v^) + b^A| (^v^) 

the triplet is given by (E[Li], 0, The density rjLi is known: 

+ {x- + {x - 

where 

(oL^ — /3^) ^ 

c(a,/3,5,A) = — 


2 K\{6y/a’^ — /?^) 


The canonical decomposition is given by. 


L( = tE[Li]-I- / / xijj — ){ds,dx) 


which has a characteristic exponent of the form 
'I'(u) = — in 

which can be obtained via: 


If ^ 

2 KxiSy'a'^ - {P + iuYY 

_\q;2 _ _|_ iu)"^ ^ 

Kx{5YcP-PA \ 


lUfJi. 


(3.27) 

(3.28) 

(3.29) 

(3.30) 

(3.31) 

(3.32) 
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GH_CF=function(u,T,r,alpha, beta, delta, nu){ 
argl=alpha* alpha-bet a*bet a 
arg2=argl-2*li*u*beta+u*u 
argm=argl-2*beta-l 

omega= - log((argl/argm)"(0.5*nu)*BesselK(delta*sqrt(argm),nu)/BesselK(delta*sqrt(argl), nu)) 
tmp=(argl/arg2)~(0.5*nu)*BesselK(delta*sqrt(arg2), nu)/BesselK(delta*sqrt(argl), nu) 
phi=exp(li*u*((r+omega)*T)+log(tmp)*T) 
return(phi) 

} 

A natural way to simulate the GH process is to directly compute, 


Lt= ^x + pYt + W{Yt). 


(3.33) 


where Yt is a random variable from a generalized inverse Gaussian (GIG) distribution. The density of the 
GIG distribution is given by, 


{x) 


2Kp{V^) 


(3.34) 


where a,b,x > 0, p G M., and Kp is a modified Bessel function of the second kind. By setting a = — fP^ 

b = S^, and p = X we obtain the density of Y. Random variables can be obtained to simulate the path of 
assets by calling: 


GH=function(alpha,beta,delta,lambda,mu, T, N){ 

h=T/N 

t=(0:T)/N 

l=rep(0, N) 

X=rep(0, N+1) 
for (i in 1:N) { 

1 [i] = rgigd, lambda, h*sqrt (alpha*alpha-beta*beta) , delta) 
X[i+1]=X[i]+mu/N+beta*l[i]+ sqrt(1[i])*rnorm(l) 

} 

return(X) 

} 


Normal Inverse Gaussian Process: By setting A = — | for the GH process, we obtain the NIG ([3]). 
The Levy triplet is inherited from before as (E[Li], 0, v^^^{dx)) where the first moment is given by. 


and the measure is defined as, 

^,Ki{a\x\)dx. 

'k\x\ 

Furthermore, the canonical decomposition is given similarly by, 

t 

Lt = t¥\Li]+ J j x{p^ — v^^^){ds,dx) 

0 R 

which has a characteristic exponent of the form, 

dt (u) = (5[-\/a2 — (/32 + iuP — \J — /3^] — iup. 


(3.35) 


(3.36) 


(3.37) 


(3.38) 
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NIG_CF=function(u,time, r,alpha,beta,delta,mu){ 

omega=delta*(sqrt(alpha*alpha-(beta+1) ~ 2 ) -sqrt(alpha*alpha-beta*beta)) 

tmp=li*u*mu*time-delta*time*(sqrt(alpha*alpha-(beta+li*u) ~ 2 ) -sqrt(alpha*alpha-beta*beta)) 
phi=exp(li*u*((r+omega)*T)+tmp) 
return(phi) 

} 

Random variables can be obtained to simulate the path of assets by calling: 

#N1G process 

NlG=function(alpha,beta,delta,mu, T, N){ 
a=l 

b=delta*sqrt(alpha*alpha-beta*beta) 

h=T/N 

t=(0:T)/N 

l=rep(0, N) 

X=rep(0, N+1) 

X[1]=0 

for (i in 1:N) { 
l[i] = lG2(a*h, b) 

X[i+1]=X[i]+mu/N+beta*delta*delta*l[i]+ delta*sqrt(1[i])*rnorm(l) 

} 

return(X) 

} 


Asset price driven by NiG process 
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Figure 8: Sample paths of an asset driven by a NIG process with o = 2, ^ = 1, (5 = 1.5, and p = 0. 


Meixner Process: Let a,(5 > 0 and —it < /3 < tt so that Li ~ MX(a,/3, 5) (|30j, |29])- Define, 



(3.39) 
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then similar to the other pure jumps models, the triplet is given by (E[Li], 0, v^^). The density rjLi is, 


VL^{x) 


(2cos(f))'' 

2anT{25) 



The canonical decomposition is given by. 


t 

Lt = tE[Li] + J j x{r]^ — iy^^){ds,dx) 

0 R 


which has a characteristic exponent of the form. 


-^{u) 


( cos(f) 

lcosh(^^); 


MX_CF=function(u,T,r,alpha, beta, delta){ 

omega=-2*delta*(logCcos(0.5*beta))-logCcos((alpha+beta)/2))) 
tmp=(cos(0.5*beta)/cosh(0.5*(alpha*u-li*beta)))"(2*T*delta) 
phi=exp(li*u*((r+omega)*T)+log(tmp)) 
return(phi) 

} 


(3.40) 


(3.41) 


(3.42) 


For simulation of the process via Brownian subordination, see |19j . Alternatively, for a simulation method 
based on acceptance-rejection sampling, see El- Random variables can be obtained to simulate the path of 
assets by calling: 

MX=function(alpha, delta, T, N) { 
h=T/N 

X=rep(0,N+l) 

for (i in 1:N){X[[i+1]]=X[[i]]+rmeixner(h,alpha,0,delta)} 
return(X) 

} 


#Meixner Increment Generator 
rmeixner = function(t,a,m,r) { 
b=0 
r=t*r 
m=t*m 
repeat { 

Q = runif(1,min=-l,max=l)/runif(1,min=-l,max=l) 

V = a/2*max(sqrt(2*r),2*r)*Q 
U = runif(1) 
if (abs(Q) < 1) { 

if (gamma(r)"2*U < abs(cgamma(r+li*V/a))“2) { 
break 

} 

} 

else { 

if (!is.na(cgamma(r+li*V/a))) 

if (max(l,2*r)*a"2*gcmima(r+l)"2*U/(2*r) < abs(cgcmima(r+li*V/a))"2 *V"2){ 
break 

} 
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} 

} 

return(V+m) 


Asset price driven by Meixner process 



Time 


Figure 9: Sample paths of an asset driven by a MX process with a = 0.5, /3 = 0, and (5 = 4. 

For a more detailed discussion on Levy processes, including their applications in finance, see m, m, m, 

HE], m, HI- 

Should there be errors in this draft, please check with published references for the correct 
formulations. 
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